Data loading:

input <- list(metpahlan_ps = "~/Documents/GitHub/oral-microbiome/data/processed/metaphlan_ps.RDS",
              motus_ps = "~/Documents/GitHub/oral-microbiome/data/processed/motus/physeq.RDS",
              pathway_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_pw_ps.RDS",
              kegg_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_kegg_ps.RDS",
              l4c_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_l4c_ps.RDS",
              go_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_go_ps.RDS",
              infogo_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_infogo_ps.RDS",
              metac_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_metac_ps.RDS",
              pfam_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_pfam_ps.RDS")

figures = "~/Documents/GitHub/oral-microbiome/figures"
dir.create(figures)

## Warning in dir.create(figures): '/Users/fconstan/Documents/GitHub/oral-
## microbiome/figures' already exists

purrr::map(input, 
           readRDS) -> objects

sample_data(objects$metpahlan_ps$physeq)$Site_Health <- paste0(sample_data(objects$metpahlan_ps$physeq)$Oral_Site,
                                                               "_",
                                                               sample_data(objects$metpahlan_ps$physeq)$Health_status)

Fig1 Identification of predominant streptococcus species at all sites

This figure is good, but with the approach we decided focusing solely on streptococcus species We should change it a bit so that it shows top 10 streptococcus species across all three sites

Identification of predominant streptococcus species at all sites

Fig1A: predominant streptococcus species across all samples in health

objects$metpahlan_ps$physeq %>%
  subset_samples(Health_status == "Healthy control") -> tmp
# subset_samples(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) -> tmp

tmp %>%
  subset_taxa(Genus %in% c("Streptococcus")) %>%
  physeq_most_abundant(physeq = .,
                       group_var = "Oral_Site",
                       ntax = 10,
                       tax_level = "Species") -> ma_saliva_tongue_strep

tmp %>%
  transform_sample_counts(function(x) x/sum(x) * 100) %>%
  phyloseq::subset_taxa(Species %in% ma_saliva_tongue_strep) %>%
  phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "Subject",#treatment
                          facet_by = c("Oral_Site", "Health_status", "Subject"),
                          ntax = 100,
                          tax_aggregate = "Species",
                          tax_add = NULL) -> p

## Warning: There are only 16 taxa, showing all

p + facet_grid(. ~ Oral_Site + Health_status, scales = "free", space = "free") + 
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 25, 50, 75, 100), 
                       labels = c(0,  0.01, 1, 10, 25,50, 75, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.1))  + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) + 
  theme(axis.text.y = element_text(angle = 0,  size = 8))-> Fig1

Fig1

Fig1 %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


Fig1 %>%
  ggsave(path = figures,
         plot = .,
         filename = "Fig1.png",
         units = "cm", width = 20, height = 10)

objects$metpahlan_ps$physeq %>%
  subset_samples(Health_status == "Healthy control") %>% 
  subset_samples(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) -> tmp

tmp %>%
  subset_taxa(Genus %in% c("Streptococcus")) %>%
  physeq_most_abundant(physeq = .,
                       group_var = "Oral_Site",
                       ntax = 10,
                       tax_level = "Species") -> ma_saliva_tongue_strep

ma_saliva_tongue_strep

##  [1] "Streptococcus_australis"     "Streptococcus_cristatus"    
##  [3] "Streptococcus_gordonii"      "Streptococcus_infantis"     
##  [5] "Streptococcus_mitis"         "Streptococcus_mutans"       
##  [7] "Streptococcus_oralis"        "Streptococcus_parasanguinis"
##  [9] "Streptococcus_salivarius"    "Streptococcus_sanguinis"    
## [11] "Streptococcus_sp_A12"        "Streptococcus_sp_F0442"     
## [13] "Streptococcus_sp_HMSC071D03"

ma_saliva_tongue_strep_sel <- c("Streptococcus_parasanguinis",
                                "Streptococcus_salivarius",
                                "Streptococcus_infantis")

objects$metpahlan_ps$physeq %>%
  subset_samples(Health_status == "Healthy control") %>% 
  subset_samples(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) -> tmp

tmp %>%
  subset_taxa(Genus %in% c("Streptococcus")) %>%
  physeq_most_abundant(physeq = .,
                       group_var = "Oral_Site",
                       ntax = 10,
                       tax_level = "Species") -> ma_saliva_plaque_strep

ma_saliva_plaque_strep

##  [1] "Streptococcus_anginosus_group"   "Streptococcus_australis"        
##  [3] "Streptococcus_cristatus"         "Streptococcus_gordonii"         
##  [5] "Streptococcus_infantis"          "Streptococcus_mitis"            
##  [7] "Streptococcus_mutans"            "Streptococcus_oralis"           
##  [9] "Streptococcus_parasanguinis"     "Streptococcus_salivarius"       
## [11] "Streptococcus_sanguinis"         "Streptococcus_sobrinus"         
## [13] "Streptococcus_sp_oral_taxon_056"

ma_saliva_plaque_strep_sel <- c("Streptococcus_oralis",
                                "Streptococcus_sanguinis",
                                "Streptococcus_gordonii")

Fig2 Overall bacterial activity of selected species

Load pathway abundance table and have a look:

objects$pathway_ps$physeq  %>%
  humann2_species_contribution(meta_data_var = c("Subject", "sample_Sample" ,"Type", "Oral_Site", "Health_status")) %>%
  separate(Feature, into = c("code", "name"), sep = ": ", remove = FALSE) %>%
  dplyr::filter(!is.na(name))  -> pwy

## Warning in speedyseq::psmelt(.): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.

## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 720 rows [1029,
## 1903, 1961, 2338, 2461, 2854, 3562, 3891, 4578, 5649, 5690, 5849, 5854, 6607,
## 6843, 7028, 7077, 7309, 7819, 8576, ...].

pwy %>%
  distinct(name, .keep_all = TRUE) %>%
  DT::datatable()

Fig2A: box plot species saliva vs. tongue

ma_saliva_tongue_strep_sel

## [1] "Streptococcus_parasanguinis" "Streptococcus_salivarius"   
## [3] "Streptococcus_infantis"

humann2_RNA_DNA_ratio_plot <- function(df,
                                       rm_un_sp = TRUE,
                                       x_plot,
                                       y_plot,
                                       color,
                                       fill,
                                       shape = NULL,
                                       filter_feature = FALSE,
                                       filter_genus = FALSE,
                                       filter_species = FALSE,
                                       only_pos = TRUE, # dplyr::filter(DNA  > 0 , RNA > 0)
                                       facet_formula = FALSE,#". ~ Health_status"
                                       export_legend = FALSE,
                                       box_width = 0.8){ #  
  
  require(tidyverse); require(ggConvexHull)
  
  if(rm_un_sp == TRUE){
    df %>%
      dplyr::filter(Species != "unclassified") -> df
  }
  
  if(filter_feature != FALSE){
    df %>%
      dplyr::filter(Feature %in% filter_feature) -> df
  }
  
  if(filter_genus != FALSE){
    df %>%
      dplyr::filter(Genus %in% filter_genus) -> df
  }
  if(filter_species != FALSE){
    df %>%
      dplyr::filter(Species %in% filter_species) -> df
  }
  
  if(only_pos == TRUE){
    df %>%
      dplyr::filter(DNA  > 0 , RNA > 0) -> df
  }
  
  df %>%
    ggplot(aes_string(x = x_plot, y = y_plot, color = color, fill = fill, shape = shape)) + 
    #geom_boxplot(outlier.colour = NA, alpha = 0.2,
                 # position = position_dodge(width=box_width)) +
    # geom_jitter(size=1, alpha=0.2) +
    # ggbeeswarm::geom_beeswarm(size=1, alpha=0.2) +
    # geom_violin(size=0.05, alpha=0.2,
    #             position = position_dodge(width=box_width)) +
    geom_hline(yintercept = 0,
               col = "red",
               linetype = "dotted",
               size = 0.5) +
    theme_light() +
    guides(fill=guide_legend(ncol=1)) +
    theme(axis.title.x = element_blank()) +
    # ylab(paste0("Relative Abundance (Top ",n," (RNA)) \n")) + #scale_y_continuous(trans='sqrt') + 
    # theme(legend.text = element_text(size= 6)) + 
    # theme(axis.text.x = element_text(size = 4)) +
    theme(legend.key.size = unit(0.2,"cm")) -> p
  
  if(facet_formula != FALSE){
    p + facet_grid(as.formula(facet_formula), drop=TRUE,
                   scale="fixed",space="free_x") + xlab(label = NULL) -> p
  }
  
  if(export_legend == TRUE){
    
    out <- list("legend" = p %>% ggpubr::get_legend() %>% ggpubr::as_ggplot(),
                "plot" = p  + ggpubr::rotate_x_text(90) + 
                  ggpubr::rotate() + theme(legend.position = "none"))
    
  }else{
    out <- p  + ggpubr::rotate_x_text(90) + 
      ggpubr::rotate() + theme(legend.position = "none")
  }
  return(out)
}

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Species", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_tongue_strep_sel,
                             facet_formula = ". ~ .",
                             export_legend = TRUE) -> plot

## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used

# plot$legend
# 
# plot$plot

ggpubr::get_palette("npg", 3)

## [1] "#E64B35FF" "#4DBBD5FF" "#00A087FF"

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

# plot$plot + geom_violin(fill = "transparent") + ggforce::geom_sina(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_violin(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#4DBBD5FF"))  + theme(legend.position = "top") -> fig2a

fig2a

fig2a %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig2a %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig2a.png",
         units = "cm", width = 20, height = 10)

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Subject", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_tongue_strep_sel,
                             facet_formula = "Oral_Site ~ Species",
                             export_legend = TRUE) -> plot

## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used

# plot$legend
# 
# plot$plot

ggpubr::get_palette("npg", 3)

## [1] "#E64B35FF" "#4DBBD5FF" "#00A087FF"

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

# plot$plot + geom_violin(fill = "transparent") + ggforce::geom_sina(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_violin(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#4DBBD5FF"))  + theme(legend.position = "top") -> fig2a

# fig2a 

fig2a %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

# fig2a %>%
#   ggsave(path = figures,
#          plot = .,
#          filename = "fig2a.png",
#          units = "cm", width = 20, height = 10)

plot$plot$data %>% 
  ggpubr::compare_means(RNA_DNA ~ Oral_Site,
                        group.by = c("Species"),
                        data = .,
                        method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05)  %>%
  DT::datatable()

Fig2C: box plot species saliva vs. plaque

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Species", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_plaque_strep_sel,
                             facet_formula = ". ~ .",
                             export_legend = TRUE) -> plot

## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

plot$plot

ggpubr::get_palette("npg", 3)

## [1] "#E64B35FF" "#4DBBD5FF" "#00A087FF"

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_violin(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#00A087FF"))  + theme(legend.position = "top") -> fig2b

fig2b

fig2b %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


fig2b %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig2b.png",
         units = "cm", width = 20, height = 10)

plot$plot$data %>% 
  ggpubr::compare_means(RNA_DNA ~ Oral_Site,
                        group.by = c("Species"),
                        data = .,
                        method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05)  %>%
  DT::datatable()

Fig 3. Pathway activity of selected species

Fig3A: significant pathway expression by candidate species saliva versus tongue in health

redo main plot + stats

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_tongue_strep_sel,
                             facet_formula = " . ~ Health_status + Species",
                             export_legend = TRUE) -> plot

## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

# plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Oral_Site,
                        group.by = c("Species","Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()

pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name",
                             y_plot = "log2(RNA_DNA)",
                             color = "Oral_Site",
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = "Streptococcus_infantis",
                             filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Health_status ",
                             export_legend = TRUE) -> plot

## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#4DBBD5FF"))  + theme(legend.position = "top") -> fig3a

fig3a

fig3a %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig3a %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig3a.png",
         units = "cm", width = 20, height = 10)

pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = "Streptococcus_parasanguinis",
                             filter_feature = diff_signif %>% filter(Species == "Streptococcus_parasanguinis") %>% pull(Feature),
                             facet_formula = " . ~ Health_status ",
                             export_legend = TRUE) -> plot

## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#4DBBD5FF")) + theme(legend.position = "top") -> fig3b

fig3b

fig3b %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


fig3b %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig3b.png",
         units = "cm", width = 20, height = 10)

pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = "Streptococcus_salivarius",
                             filter_feature = diff_signif %>% filter(Species == "Streptococcus_salivarius") %>% pull(Feature),
                             facet_formula = " . ~ Health_status ",
                             export_legend = TRUE) -> plot

## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp


ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#4DBBD5FF"))  + theme(legend.position = "top")  -> fig3c

fig3c

fig3c %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig3c %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig3c.png",
         units = "cm", width = 20, height = 10)

Fig3B: significant pathway expression by candidate species saliva versus plaque in health

redo main plot + stats

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_plaque_strep_sel,
                             facet_formula = " . ~ Health_status + Species",
                             export_legend = TRUE) -> plot

## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

# plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Oral_Site,
                        group.by = c("Species","Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()

pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name",
                             y_plot = "log2(RNA_DNA)",
                             color = "Oral_Site",
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = "Streptococcus_oralis",
                             filter_feature = diff_signif %>% filter(Species == "Streptococcus_oralis") %>% pull(Feature),
                             facet_formula = " . ~ Health_status ",
                             export_legend = TRUE) -> plot

## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#00A087FF"))  + theme(legend.position = "top") -> fig4a

fig4a

fig4a %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig4a %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig4a.png",
         units = "cm", width = 20, height = 10)

pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = "Streptococcus_sanguinis",
                             filter_feature = diff_signif %>% filter(Species == "Streptococcus_sanguinis") %>% pull(Feature),
                             facet_formula = " . ~ Health_status ",
                             export_legend = TRUE) -> plot


plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#00A087FF"))  + theme(legend.position = "top") -> fig4b

fig4b

fig4b %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig4b %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig4b.png",
         units = "cm", width = 20, height = 10)

# pwy %>%
#   dplyr::filter(Health_status == "Healthy control") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Oral_Site", 
#                        fill = "Oral_Site",
#                        shape = NULL,
#                        filter_species = "Streptococcus_gordonii",
#                        filter_feature = diff_signif %>% filter(Species == "Streptococcus_gordonii") %>% pull(Feature),
#                        facet_formula = " . ~ Health_status ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Fig 5. Overall bacterial activity of selected species in health versus diseases

Fig5A: box plot species in health and disease at each site

strep_sel <- c(ma_saliva_plaque_strep_sel, ma_saliva_tongue_strep_sel)

pwy %>%
  # dplyr::filter(Health_status == "Healthy control",
  # Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Species", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = strep_sel,
                             facet_formula = ". ~ Oral_Site",
                             export_legend = TRUE) -> plot

## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

plot$plot + ggpubr::rotate()

plot$plot$data %>% 
  ggpubr::compare_means(RNA_DNA ~ Health_status,
                        group.by = c("Oral_Site", "Species"),
                        data = .,
                        method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05)  %>%
  DT::datatable()

ggpubr::get_palette("jco", 3)

## [1] "#0073C2FF" "#EFC000FF" "#868686FF"

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_violin(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#EFC000FF", "#868686FF"))  + theme(legend.position = "top") -> fig5

fig5

fig5 %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


fig5 %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig5.png",
         units = "cm", width = 20, height = 10)

Fig 6 Pathway activity of selected species at each site in health versus disease

Fig6A-F: significant pathway expression by candidate species at each site in health and disease

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  # dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             # filter_species = "Streptococcus_infantis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species + Oral_Site ",
                             export_legend = TRUE) -> plot

plot$legend

# plot$plot

A - SALIVA

strep_sel <- c(ma_saliva_plaque_strep_sel, ma_saliva_tongue_strep_sel)
strep_sel

## [1] "Streptococcus_oralis"        "Streptococcus_sanguinis"    
## [3] "Streptococcus_gordonii"      "Streptococcus_parasanguinis"
## [5] "Streptococcus_salivarius"    "Streptococcus_infantis"

Streptococcus_infantis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_infantis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot

# plot$legend
# plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") -> diff

diff %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()

# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_infantis") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_parasanguinis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_parasanguinis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot

# plot$legend
# plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()

pwy %>%
  dplyr::filter(Species %in% "Streptococcus_parasanguinis") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             # filter_species = "Streptococcus_infantis",
                             filter_feature = diff_signif %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot

## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#EFC000FF", "#868686FF"))  + theme(legend.position = "top") -> fig6a

fig6a  %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig6a %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig6a.png",
         units = "cm", width = 20, height = 10)

Streptococcus_salivarius

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_salivarius",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot

# plot$legend
# plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()

pwy %>%
  dplyr::filter(Species %in% "Streptococcus_salivarius") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             # filter_species = "Streptococcus_infantis",
                             filter_feature = diff_signif %>% pull(Feature),
                             facet_formula = "  ~ Species ",
                             export_legend = TRUE) -> plot

## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#EFC000FF", "#868686FF"))  + theme(legend.position = "top") -> fig6b


fig6b  %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


fig6b

fig6b %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig6b.png",
         units = "cm", width = 20, height = 10)

Streptococcus_gordonii

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_gordonii",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot

plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()

# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_gordonii") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_sanguinis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_sanguinis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot

plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()

# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_sanguinis") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_oralis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_oralis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot

plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()

pwy %>%
  dplyr::filter(Species %in% "Streptococcus_oralis") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             # filter_species = "Streptococcus_infantis",
                             filter_feature = diff_signif %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot

## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used

plot$legend

plot$plot

B - TONGUE_BIOFILM

strep_sel <- c(ma_saliva_tongue_strep_sel)
strep_sel

## [1] "Streptococcus_parasanguinis" "Streptococcus_salivarius"   
## [3] "Streptococcus_infantis"